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A differential analysis of the power flow equations 
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Abstract — The AC power flow equations are fundamental in 
all aspects of power systems planning and operations. They are 
routinely solved using Newton-Raphson like methods. However, 
there is little theoretical understanding of when these algorithms 
are guaranteed to find a solution of the power flow equations or 
how long they may take to converge. Further, it is known that in 
general these equations have multiple solutions and can exhibit 
chaotic behavior. In this paper, we show that the power flow 
equations can be solved efficiently provided that the solution 
Ues in a certain set. We introduce a family of convex domains, 
characterized by Linear Matrix Inequalities, in the space of 
voltages such that there is at most one power flow solution in 
each of these domains. Further, if a solution exists in one of these 
domains, it can be found efficiently, and if one does not exist, a 
certificate of non-existence can also be obtained efficiently. The 
approach is based on the theory of monotone operators and 
related algorithms for solving variational inequalities involving 
monotone operators. We validate our approach on IEEE test 
networks and show that practical power flow solutions He within 
an appropriately chosen convex domain. 

I. Introduction 

Power systems are experiencing revolutionary changes 
due to various factors, including: Integration of renewable 
generation, distributed generation, smart metering, direct 
or price-based load-control capabilities. While potentially 
contributing to the long-term sustainability of the power 
grid, these developments also pose significant operational 
challenges by making the power system inherently stochas¬ 
tic and inhomogeneous. As these changes become more 
widespread, the system operators will no longer have the 
luxury of large positive and negative reserves. Moreover, 
operating the future power grid will require developing new 
computational tools that can assess the system state and 
its operational margins more accurately and efficiently than 
current approaches. Specifically, these new techniques need 
to go beyond linearized methods of analysis and ensure 
that the power system is safe even in the presence of large 
disturbances and uncertainty. In this paper, we focus on the 
fundamental equations of the power system - the Power 
Flow (PF) equations. The equations constitute a system 
of nonlinear equations and are known to exhibit complex 
and chaotic behavior [1][2]. In the past, when the power 
systems were operated well within their security margins, the 
PF equations were solved without difficulty using standard 
techniques like Newton-Raphson and its variants. However, 
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changes in the power system mentioned above mean that 
this may no longer be possible. Thus, Newton-Raphson tech¬ 
niques which rely on good initialization may fail to converge. 
In such a situation, it becomes difficult to assess whether 
the system is actually operationally unsafe or the Newton- 
Raphson method failed because of numerical difficulties or 
bad initialization. In this paper, we propose an approach to 
remedy this problem. Our approach is based on the theory of 
monotone operators. Just as a convex optimization problem 
can be solved efficiently, one can find zeros of a monotone 
operator efficiently. (In fact, the gradient of a convex function 
is a monotone operator.) Thus, if we can show that the 
nonlinear PF equations can be described by a monotone 
operator, then they can be solved efficiently. It turns out 
that the PF operator is not globally monotone, however it 
is monotone over restricted domains. 

Our main technical result is a description of the do¬ 
main over which the power flow operator (whose zeros are 
solutions to the PF equations), is monotone. In fact, we 
introduce a family of monotonicity domains, parameterized 
by a matrix-valued parameter. Each monotonicity domain is 
characterized by a Linear Matrix Inequality (LMI) in the 
real and imaginary components of the voltage phasor at each 
bus. Within each of these monotonicity domains, there can 
be at most one solution of the PF equations. Further, the 
solution of the PF equations within each domain can be 
reduced to the solution of a monotone variational inequality, 
for which there exist efficient algorithms. The algorithms run 
in polynomial time and terminate either by finding a solution 
to the PF equations or a certificate that no solution exists to 
the PF equations within that domain. 

The choice of monotonicity domain is critical, since dif¬ 
ferent domains cover different parts of the voltage space. 
Ideally, one would like to find a monotonicity domain that 
covers all the solutions of interest, i.e. solutions that satisfy 
typical operational constraints on voltages and phases. In 
general, this is a hard problem. We deal with this problem by 
proposing a sufficient condition to ensure that voltages satis¬ 
fying some bounds are constrained within the monotonicity 
domain. Finding the largest bound can then be cast as a 
quasi-convex optimization problem. We use this technique to 
pick the monotonicity domain. Numerical tests show that this 
approach is able to generate a sufficiently large monotonicity 
domain for several test networks. 

The rest of this paper is organized as follows. Section |I^ 
covers relevant background on power systems and monotone 
operators. The main technical results are presented in Section 
m In Section |IV] we discuss how our approach compares 
to related work on conditions for existence and uniqueness 


of PF solutions. In Section |Vj we present numerical results 
illustrating our approach on some IEEE benchmark networks. 

II. Modeling Power Systems 

A. Notation 

R is the set of real numbers, C the set of complex 
numbers. R”,C” denote the corresponding Euclidean space 
in n dimensions. 

Given a set C C R", Int (C) denotes the interior of the 
set. Given a complex number a; S C, Re (x) denotes its real 
part and Im (x) its imaginary part, x* denotes its complex 
conjugate. If X G is a square matrix with complex 

entires, X* denotes the conjugate transpose. ||a;|| refers to 
the Euclidean norm of a vector x G R" or a; S C" and 
{x, y) to the standard Euclidean dot product. Given an vector 
x G R", di (a;) denotes the nxn diagonal matrix with {i, i)- 
th entry equal to Xi. denotes the space of kxk symmetric 
matrices. 

Given a differentiable function / : R^ >->• R^, V/ denotes 
the Jacobian of f, a k x k matrix with the z-th row being 
the gradient of the z-th component of /. 

Eor M G R"^”, Sy(M) = M+M^_ j the indicator 
function; 

1 if p = True 
0 if p = Ealse 

B. Power System Model 

The transmission network is modeled as a graph iV^E) 
where V is the set of nodes and E is the set of edges. In 
power systems parlance, the nodes are called buses and the 
edges are called lines (transmission lines). We shall use these 
terms interchangeably in this manuscript. Nodes are denoted 
by indices z = 0,1,...,rz and edges by ordered pairs of 
nodes (z,j). We pick an arbitrary orientation for each edge, 
so that for an edge between z and j, only one of (z,j) and 
(j,z) is in E. If there is an edge between buses z and j, we 
write z ~ j, j ^ i. 

Edges correspond to transmission lines. The transmission 
network is characterized by its complex admittance matrix 
Y G C"^”. Y is symmetric but not necessarily Hermitian. 
Let G = Re (y), S = Im (F). 

Let Vi be the voltage phasor. Pi and Qi denote active and 
reactive injection at the bus z respectively. V is the vector of 
voltage phasors at all buses. Buses are of three types; 

• (P,V) buses where active power injection and voltage 
magnitude are fixed, while voltage phase and reactive 
power are variables. The set of (P,V) buses is denoted 
by pv. The voltage magnitude setpoint at bus z G pv is 
denoted by Vi. 

• (P,Q) buses where active and reactive power injections 
are fixed, while voltage phase and magnitude are vari¬ 
ables. The set of (P,Q) buses is denoted by pq. 

• Slack bus , a reference bus at which the voltage mag¬ 
nitude and phase are fixed, and the active and reactive 
power injections are free variables. The slack bus is 
denoted by S and its voltage phasor by Vq- We choose 
bus 0 as the slack bus as a convention. 


C. Background 

1) Power Flow Equations: The PE equations model the 
flow of power over the network. They are a set of coupled 
nonlinear equations that follow from Kirchoff’s laws applied 
to the AC power network. Circuit elements in the standard 
power systems models are all linear, if one ignores discrete 
elements like tap-changing transformers. Even though Ohm’s 
laws and Kirchhoff’s laws are linear in voltages and currents, 
power is a product of a voltage and a current and hence 
quadratic. 

PE equations are static and as such the equations model the 
regime when the network is balanced, that is, the net sum of 
power consumptions, injections and power dissipated is zero. 
This relies on the assumption that at the time-scale where 
the PE equations are solved (every few minutes), the system 
is in a quasi-steady state, i.e. the dynamic disturbances 
have been resolved through actions of the automatic control 
(voltage regulators, power system stabilizers and primary and 
secondary frequency control systems). The complete set of 
PE equations over the graph (V,£) is formally stated as 

Pi = Re {Vi{YV)*) z G pq U pv 
Qi = fm(yi{YV)*) z G pq 

\Vi\=Vi i G pv 
Ks = Vo 

It will also be convenient, for what follows, to utilize the 
Cartesian parametrization of voltages; 

V = V,^+iV^^,iGpqUpv. 

fv^\ 

Let V‘^ = \yy] denote the stacked vector of real and 
imaginary voltage components, so that = 

Definition 1. Define the power flow operator as F : R^" i— 

R2n 

[F(y‘^)V = Gu({ytf + {vn'') 

jn^i 

- ^ G,, {v-v- + y/y/) - p,, z G pv u pq 

jr^i 

(la) 

[p(y^)]„+, = p..((ynV(y/f) 

+ Y^B,,{vfvf + vyvy) 

jr^i 

+ E {y^y^ - y^y^) -Q^p^ pq 

jr^i 

(lb) 

[F iV^)U, = (Vff + - vl z G pv (Ic) 

Then the PE equations can be written as 

p(y'=) = 0. 






D. Monotone Operators 

We now review briefly the theory of monotone operators, 
as is relevant to the approach developed in this paper. For 
details and proofs of the results quoted in this Section, we 
refer the reader to the recent survey [3]. A function H : 

I—is said to be a monotone operator over a convex 
domain C if 

{H{x)- H{y),x-y) >0 'ix,yeC 

A monotone operator is a generalization of a monotonically 
increasing function (indeed, if A: = 1, the above condition 
is equivalent to monotone increase: x > y H {x) > 

H (y)). H is said to be strongly monotone with modulus m 
(or simply strongly monotone) if 

{H{x)-H{y),x-y)>'^\\x-yf \/x,yGC 

for some m > 0. A common example of a monotone operator 
is the gradient of a differentiable convex function. 

Definition 2 (Monotone Variational Inequality). Let C C 
be a convex set and H he a monotone operator over C. The 
variational inequality (VI) problem associated with H and C 
is: 

Find x & C such that {H {x),y — x) >0 yy G C (2) 

The following result shows that monotone variational 
inequalities with compact domains always have a solution 
and can be solved efficiently. 

Theorem II.l. If H is strongly monotone operator over 
a compact domain C, then has a unique solution x*. 
Further, an approximate solution x,, G C satisfying 

||a;£ - a:|| < e (3) 

can be found using at most O (log ( 7 )) evaluations of H 
and projections onto C. 

Remark 1. In this manuscript, we are interested in finding 
zeros of the PF operators introduced above. We can use 
monotone operator theory for this as follows: Suppose H 
satisfies the hypotheses of theorem [lI.l| If there exists a point 
X* G C with H {x*) = 0, then this is the unique solution 
of the variational inequality. Conversely, if the variational 
inequality has a solution with H (x*) 7 ^ 0, then have a 
certificate that there is no solution of H (x) = 0 with x G C. 

The next result provides a simple characterization of 
monotonicity for differentiable operators: 

Theorem II.2. Suppose that H is differentiable. Then H is 
strongly monotone with modulus m > 0 over C if and only 

if 

VH{x)AVH{xy >ml Vx S C 
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monotone, the PF solutions can be found efficiently. Since 
PF equations can have multiple isolated solutions, it is not 
possible that the PF operator is globally monotone because 
this would imply a unique solution to the PF equations. 
Thus, we focus on characterizing domains over which the 
PF operator (or a scaled version of it) is monotone. Proofs 
of all theorems in this section are deferred to the appendix 
Section |Vl] 


A. Characterization of Domain of Monotonicity of the Power 
Flow Operator 

We now derive our main results on the monotonicity of 
the PF operator Q. In order to state the result succinctly, 
we will need to define some matrices that are functions of 
the network topology, locations of (P,V) buses and of the 
network impedance matrices. 

Definition 3. Define the following row vectors: 

G* = (Gil, Gi 2 , . . . , Gin) 

= {Bil,Bi2, . . . , Bin) 

Gpq = (G^I (1 G pq),..., Gi„I (n G pq)) 

■Spq = (Sill (1 G pq),..., Binl (n G pq)) 


Let Ci denote the Lth column of the n x n identity matrix, 
so Ci has zeros everywhere except the i-th entry. 

Definition 4. Define the following matrices: 


M, 


N, 


Mo 

No 


' f di (G*) + e,G* di (S*) - e,B^ \ 

^ \-^i{BU)-GBl^ di(Gg-e,G*J 

/ di (G*) + e,G* di (S*) - CiSA 
^ V“di (S(,q) + 2eie,'^ di (G(,q) j 


I -di (S*) + CiB^ di (G*) + CiG* \ 

-di (G(,q) + e.G(,q -di (B^q) - c.B^) 

-di (B*) + CiB^ di (G*) + e,G* \ 
-di (G(,q) -di (S^q) + 2e,e,^y 

( di (GO) di {B°) \ 

l-di(S°q) di(GOq); 

( -di (5°) di (GO) \ 
l-di(G°q) -di(SOq)j 


i G pq 

i G pv 
(4a) 

i G pq 

i G pv 

(4b) 

(4c) 

(4d) 


We now state our main technical result, which shows that 
the PF equations can be solved by solving a monotone varia¬ 
tional inequality (for which there exist efficient algorithms). 

Theorem III.l. Let W G K 2 nx 2 n invertible matrix 

and m > 0. There is at most one solution of the PF equations 
F (L°) = 0 over the domain: 


In this Section, we study the monotonicity of the PF 
operator F Q. As described in Section |II-D| zeros of F 
(solutions to the PF equations) can be found efficiently if F 
is monotone. Thus, if we can prove that the PF operator is 


C{W) = 



: ^ Sy (IL {M,Vf + N,Vf)) h ml \ 

i^O ) 

(5) 









Define Fw (V'^) = WF(V‘^). Fw is strongly monotone 
with modulus m over the set C {W). Let V‘^* be the unique 
solution to the monotone variational inequality: 

&C [W] (6) 

{Fw (F^), a - > 0 ya-.a&CiW) (7) 

If F = 0, V^* is the unique solution to the PF 

equations in C {W). Otherwise, there are no solutions to 
F{V^) = Q,V<^ &C{W). 


Remark 2. We implicitly assume existence of a solution 
to the variational inequality above. This is guaranteed (by 
theorem II. li if C (W) is compact. Magnitudes of voltage 
phasors in practical power systems are always bounded, and 
one can simply pick a large positive number 6 > 0 such that 
all power flow solutions satisfy ||F'^|| < b and define C {W) 
with this additional constraint, ensuring it is a compact set. 

1) Illustration on 2-bus network: We consider a 2-bus 
network with admittance matrix 


/.05-jl.ll -.05-fjl.ll\ 

l^-.05-f jl.ll .05-jl.ll ) 


We flx = /,I/o = 1 +j0. Let Vi - Vb = By 

varying W, we can And two disjoint monotonicity domains 
for this system: 


Wi = 


f-6.75 

I 0 



f-6.75 

V 0 



It can be verified that these two monotonicity domains 
cover the space of {V^, V^) and there is at most one solution 
in each domain. In fact, in this system, as long as there is 
a solution, there exist two (one in each domain). As the 
loading increases, the two solutions move closer to each 
other and eventually disappear (at the point of the saddle- 
node bifurcation). Note that the boundary between the two 
solutions is observed exactly at Vi — Vb = Note that 
letting Vi = V exp (j0), or equivalently 

V cos (9) > i 

which is a well-known classical criterion for voltage stability 
in the two bus example. 


B. Selection of the Monotonicity Domain 


Theorem III.l shows that PF solutions in C (W) can be 
found efficiently. However, the matrix inequality characteriz¬ 
ing C {W) is not intuitive and does not have a simple inter¬ 
pretation. Further, the choice of W would affect the domain 
of monotonicity. Thus, it is important to pick W so that 
C {W) contains the “operationally relevant” PF solutions. In 
this Section, we describe a technique to achieve this goal. 

We consider the following class of operational constraints: 


\Vi - Vi\ < 6i 


where u is a “nominal voltage profile. Intuitively, one can 
think of this as the “center” of the monotonicity domain. 


Stated in terms of = 


yx 

yv 


, the constraint becomes 


- uf p -f |y/ - ufp < 6i i = l,...,n 

Cop {5) = 

{ (yy) ■ * = 


( 8 ) 


We now state a theorem that gives a sufficient condition 
to guarantee that Cop (<)) C C (W). 


Theorem III.2. Let 

' 1-^2 -1 -1 1 - v / 2 ' 

1 v/2-1 I-V 2 -1 

Suppose 3W, Xi ,..., Xn G 5^" satisfying 


K = 


Y, Sy {W {M,vf + )) hml + Yx^ 

i—0 2—1 

X, h 4 i-KuM, - K2iNi) ,; G {1, ..., 4} 
^ 4 {KuM, - K^iNi) ,; G {1,... , 4} 


(9) 

( 10 a) 

(10b) 

( 10 c) 


Then Cop (^) C C (W). Further, this is a convex feasibility 
problem in W, Xi,, Xn. 

Theorem |III.2| gives us a principled way to choose W so 
that Cop C C iW). However, the system of constraints ( [T0| ) 
may be infeasible. In this situation, we can try to see if a 
scaled version of the constraint set C (b) lies within C {W). 
Specifically, we look for the largest possible p > 0 such that 
C {p5) C C (W). 

We then formulate the following problem: 

Maximize p (Ha) 

Subject to 

n n 

^ Sy (IF ^ ml(11b) 

2—0 2—1 

X, ^ p4 {-KiiM, - KM,I G {1,... ,4} (11c) 

X,hpSk{KuM,-K2iNi),lG{l,...,A} (lid) 


lemma 1. (|1 l| i is a quasi-convex optimization problem and 
can be solved efficiently. Further, if the matrix 

n 

i=0 

is full-rank, the problem is always feasible and the optimal 
solution satisfies Cop {pS) C C (IF). 

Remark 3. The condition that X]"=o is non¬ 

singular means that the Jacobian of the power flow equations 
at the nominal voltage profile v is non-singular. This is 
generically true, except for very special choices of the v 
and the network admittance matrix. Further, power system 
voltage stability criteria are often formulated in terms of the 
“distance to singularity” of the Jacobian matrix, so enforcing 
a non-singular Jacobian is a reasonable restriction. 








IV. Related Work 

Several papers have studied conditions for existence of 
solutions to the PF Equations [4][5][6]. In [7], the authors 
propose a sufficient condition for the insolvability of the PF 
equations based on a convex relaxation. Our approach differs 
from these in the following important ways: 

• To solve the PF equations in C {W), we provide neces¬ 
sary and sufficient conditions, that is, our approach finds 
a solution in C {W) if there exists one, and a certificate 
of non-existence if there is no solution. 

• Our approach is algorithmic, that is, we provide an 
algorithm (based on solving a monotone variational 
inequality) that is guaranteed to find the solution ef¬ 
ficiently, i.e. in polynomial time. 

• If there are additional operational constraints 

yy) e S, where S' is a convex set and S C C [W) 
, e.g. corresponding to line protection [8] and/or line 
flow limits, we can additionally answer the question of 
whether there exists a PF solution with y^) G S. 
This is an important contribution, since most of the time 
system operators are interested in finding PF solutions 
that additionally satisfy operational constraints. 

• We can And multiple PF solutions by choosing differ¬ 
ent nominal voltage profiles v. This is relevant for a 
number of important power system applications, such 
as assessing distance to voltage collapse and transient 
stability (computation of the so-called controlling un¬ 
stable equilibrium point). 

In [4], the authors also propose an algorithm based on a 
contraction mapping. However, the algorithm only works in 
a small ball around the origin in the {P, Q) space. This was 
extended to other kind of sets in [9]. On the other hand, our 
results are stated in terms of a convex constraint in (V®, V^) 
space. Understanding the set of (P, Q) for which the solution 
(y^yy) g c {W), and conversely the set of (V^yy) for 
which {P, Q) lies in a certain set, is still an open problem, 
even for the special case where all buses are (P,V) buses. 
This setting was studied and the results were extended 
and connected with results on synchronization in coupled 
oscillators in a series of recent papers [10][11][12]. These 
authors provided distinct sufficient and necessary conditions 
on the injections for the existence of power flow solutions 
with phase differences satisfying certain bounds. 

In recent work [13], we have shown that for the special 
case of lossless networks (and networks with constant ratio 
of inductance to resistance) the PF equations can be solved 
analyzing a convex optimization problem. Specifically, we 
minimize the so-called energy function over a restricted 
domain (the convexity domain of the energy function). 
The optimization problem in [13] was formulated in polar 
coordinates, i.e. in terms of the voltage magnitude and 
phase at each bus. In that setting, we get a single domain 
characterized by a nonlinear but convex matrix inequality, 
and any solution within this domain can be found efficiently. 
In fact, if one writes the PF equations in polar coordinates, 
one can show that the monotonicity domain coincides pre¬ 


cisely with the convexity domain of the energy function. 
However, this polar-coordinate monotonicity domain of [13] 
is not equivalent to monotonicity domain in the Cartesian 
coordinates discussed in this manuscript. In this context, 
work reported in the present manuscript was inspired by 
our earlier attempt to extend the polar-coordinate based 
approach to lossy networks. Even though we were not able 
to And a simple convex characterization of the monotonicity 
domain in the polar coordinates, it lead us to discover that 
in the Cartesian coordinates the monotonicity domain can 
be described using LMIs. Numerical tests show that the 
domain of convexity of the energy function is a superset 
of the monotonicity domain computed by for lossless 
networks. Furthermore, in that setting, we were able to prove 
a stronger result for tree networks showing that there exists 
a PF solution if and only if there exists one in the convexity 
domain. The condition characterizing the convexity domain 
also had a simple physical interpretation - it simply requires 
voltage magnitudes and phases at the neighboring buses to 
be “close” to each other. 

However, the approach reported in the present manuscript 
is more general. It allows one to tune the monotonicity 
domain to compute multiple solutions. Furthermore, the 
monotonicity domain is expressed as an LMI here and can 
be handled using off-the-shelf software, while the nonlin¬ 
ear matrix inequality from [13] requires specially designed 
solvers. Understanding the relationship between the approach 
of this manuscript and the one from [13], in particular, the 
relationship between cartesian and polar parameterizations, 
is an important direction for future work. 

V. Numerical Illustrations 

This Section presents numerical experiments illustrating 
theoretical results presented above. We start by illustrating 
monotonicity domain on a 3 bus system example, and 
then discuss numerical experiments performed on the case9 
and casel4 systems available with the MATPOWER [14] 
software. For solving the convex optimization problem 
we use the parser-solver CVX [15][16]. For solving the 
variational inequality, we implement our own solver based 
on the extra-gradient method described in [3]. The imple¬ 
mentations used here are not optimized and currently do not 
scale to larger systems. Algorithmic developments that would 
enable applications of these ideas to larger network will be 
an important direction for future work. 

A. Illustration of Monotonicity domain for 3-bus network 

We consider a 3 bus network with bus 0 being the slack 
bus, bus 1 a (P,V) bus and bus 2 a (P,Q) bus. The voltage 
phasor at the slack bus is taken to be 1 -f jO and at the 
(P,V) bus to be exp(j0). The phasor at the (P,Q) bus is 
parameterized as 1 -f jO -f where represent 

the deviations from the reference voltage. The variables to 
be solved for in the PF equations are V^yy,9. We pick 
a nominal voltage profile by setting all voltages equal to 
1 -f jO. We use (ED to And a monotonicity domain W. We 


then plot the monotonicity domain in space for 

multiple values of 9 in Fig. Q- 

Fig. Q shows that for small values of 6, the monotonicity 
domain covers a fairly large space in V^. As 6 increases, 
the domain shrinks and ultimately, before |0| hits it 
becomes empty. This is consistent with the idea that all 
voltages “close” to the nominal voltage profile are contained 
within the monotonicity domain. Fig. Q also shows that 
the domain actually includes a fairly large region around the 
nominal voltage profile. 



(a) 61 = 0 (b) e = .23tt 



(e) 9 = —.457r 

Fig. 1: Monotonicity Domain for a 3 bus system 


B. Computation of Voltage Phasor Bounds 


We use the methodology of Section III-B to find the 


maximum possible deviation from a nominal voltage profile 
contained within a monotonicity domain. We choose Vi = 
Vo, so that every bus has a voltage phasor equal to the 
reference at the slack bus. For a system with only (P,Q) buses, 
this is in fact a solution to the PF equations for the case of 
no, 0, injections. In practical power systems, large deviations 
from the slack bus voltage are relatively rare. Therefore, it 
is reasonable to seek a monotonicity domain that includes 
all solutions in some neighborhood of this special (nominal) 
one. 

We choose a uniform bound 6^ = 1 and solve Q to com¬ 
pute W. We report the bounds for our test networks in the 
Table These bounds easily suffice to cover typical voltages 
observed in practical power systems operations. However, 
these domains may be insufficient to describe solutions in 
the case of extreme injections used in planning studies and 
contingency analysis. Our numerical experience suggests that 
these bounds are are conservative. More specifically, there 
are many that do not satisfy the bound given by 

p, but still lie within C {W). The results of the next section 
establish this numerically. 


System 

Limit on |Vi — Vb| (p.u system) 

Case 9 

.2 

Case 14 

.25 


TABLE I; Bounds on guaranteeing V‘^ £ C (W) 


C. Scaling Loads 

In this Section, we use the monotonicity domain obtained 
by solving but instead of relying on the computed 

bound p, we test numerically whether existence of a solution 
of the PF equations implies existence of a solution in the 
monotonicity domain. Of course, determining existence of 
a solution to the PF equations is a hard problem, so we 
instead rely on sufficient conditions for insolvability devel¬ 
oped in [7]. We use the implementation available with the 
MATPOWER package [17]. 

Our experimental approach is as follows: We take the 
injections given as part of the Matpower test case, and scale 
each injection by a complex scalar: Si = aSi. By scaling 
the magnitude of injections, we reach a point where there 
are no solutions to the power flow equations (maximum 
loadability). By scaling by a complex factor, we effectively 
achieve different scaling for the real and reactive parts of the 
injections. We pick the magnitude of a uniformly distributed 
between 1 and a maximum value and the phase uniformly 
distributed between 7r/6 and 7r/3. (We observe that once 
we scale beyond a certain threshold there are very few 
cases where the PF equations have a solution.) These choices 
create a space of injections that include at least part of 
the solvability boundary (injections beyond which the PF 
equations have no solution). For a uniform grid of the a 
space described above with 100 points, we compute the 
number of points at which there are no PF solutions (that 
is in this case we have a certificate of infeasibility based on 
[7]) and the number of points at which there are no solutions 
in C {W). The results are shown in Table 0 - 


System 

C[W) 

-C{W) 

NoSol 

Case 9 

38 

6 

56 

Case 14 

40 

49 

11 


TABLE II: Existence of Solutions in the Monotonicity Do¬ 
main: The first column denotes the number of instances 
(out of 100) such that a PF solution within C {W) was 
found. The second denotes the number of instances for which 
a PF solution within C (W) was not found, and there is 
no certificate of insolvability of the power flow equations 
based on the sufficient condition from [7]. The last column 
denotes the number of instance where we have a certificate 
of infeasibility. 


The results show that for the 9 bus network, most cases 
where a solution exists lie within the monotonicity domain. 
However, for the 14 bus system, there exists a significant 
number of cases for which the solution does not lie in the 
monotonicity domain, if it does exist. It turns out that there 
are larger monotonicity domains that contain the original one 



















(we checked this post-hoc by solving a feasibility problem 
to find a W such that C {W) includes all the solutions found 
using Newton-Raphson). Thus, it seems like the choice of W 
based on GD is not optimal and larger monotonicity domains 
could be found. How one might do this in a principled way 
constitutes a comprehensive direction for future work. 

D. Comparison to Energy Function Approach 

In this Section, we consider a lossless modification of the 9 
bus and 14 bus networks and perform the same experiment as 
described in the previous Section. This allows us to compare 
the approach developed in this manuscript with the energy 
function based approach described in [13]. 

The results (table |In| ) show that for these networks, exis¬ 
tence of a solution seems to imply existence within Cef but 
not within Cw- This suggests that Cef is a superset of 
and that potentially choosing W based on ( [TT] ) does not cover 
all the solutions. Closing this gap and designing approaches 
to ensure that Cw would contain all the solutions in Cef is 
a direction for future work. We anticipate that development 
of this approach will likely suggest improvements applicable 
to lossy networks as well. 


System 

C[W) 

Cef 

-C(W) 

—Cef 

NoSol 

Case 9 LL 

50 

56 

6 

0 

44 

Case 14 LL 

62 

92 

30 

0 

8 


TABLE III: Monotonicity vs Energy Eunction: The first, 
third and final columns have the same meaning as table 0 - 
The second column is the number of instances where a PE 
solution within Cef was found, and the final column the 
number of instances where a PE solution within Cef was 
not found and there is no certificate of infeasibility. 


VI. Conclusions and Euture Work 

We have presented a novel approach to solving the PE 
equations based on monotone operator theory. Our main 
technical contribution is a characterization of the family of 
domains over which the PE operator is monotone. Within 
any of these domains, there can be at most one power 
flow solution. The approach leads to efficient algorithms for 
determining existence and actually finding solutions to the 
PE equations within the monotonicity domains. In spite of 
the progress reported, there still remain many unanswered 
questions. Eirst, the numerical experiments show that our 
approach to selecting W based on solving GD is potentially 
conservative and misses some solutions that could be cap¬ 
tured by other choices of W. Eiguring out the relationship 
between the monotone-operator approach developed h‘ere 
and the energy function approach of [13] is another important 
direction for future work. Einally, developing algorithms that 
scale well and allow these techniques to be applied to larger 
systems is another item in our path forward. 
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Appendix 

A. Proof of theorem m 

Proof. The result is a direct consequence of the strong 
monotonicity of Fw over C {W). By lemma the Jacobian 












of Fw is given by 


VFw W) = Y.^ 

z=0 


where /3 = 


yy)' 


F]y is strongly montone with modulus 


if (theorem |II.2|) 


Sy (VFw m = ^ Sy + N^V^y)) F ml 

i=0 


Thus, Fyy is strongly monotone over the domain f3 G C (W). 
The remaining results follow from theorem |II.1| □ 

B. Proof of theorem |///■2| 

Proof Let Qq = YJi=o Sy {W {M,vf + N.yy)). Then, 

n 

Y,Sy{W{M,Vf + mvy)) = 

n 

Qo + Sy {W (M, {Vf - <) + TV, (Vf - vf))) (12) 

2=0 


e Cop{s) = ^\vf-vf\^ + \vy-vf\^<s. 



This can be seen simply by plotting the points 
Si{Kii+jKi 2 ) on the complex plane and observing 
that their convex hull contains the disc of radius Si centered 
at the origin. Thus, for any 1^° G Cop (5), 


-{M,ivf-vt) + N,{vy-vy))G 

<5,Conv ({-iVyM, - K 2 iN,,KnM, - 
From ( [T^ , we now have 

n 

Y, Sy (W^ {M,Vf + N^Vf)) hQo-Y ^ 

2=0 k^£ 

By ( |llb| l, we have G Cop (S) G C (W), so that 

Cop(S)cC(lV). □ 

lemma 2. The Jacobian of F {V^, Vy) is given by: 

n 

\/F{V^,Vy)=Y^i^f+ (13) 

2=0 

Proof. The entries of the Jacobian are given by; 

vF{v^,vy)=(^^ 

with each block being n x n. 

Su = 2G,,Vf + Y - E 

i~* j~2 

S^J = - G,,Vr 

Tu = 2G„y/ - 5] - Y 

jr^i j^i 

7), = B,,Vf - G.,vy 


The entries of O, L depend on whether i is a (P,Q) bus or 
(P,V) bus. For i G pq, we have 

0„ = 2B,,Vi^ + Y - Y G^Jyf 

j^i jr^i 

o,, = B,,Vf + 

L„ = 2BuVy + E - Y G^yf 

j^i jr^i 

= B^jVf + G,jVf 
For i G pv, we have 

0„ = O,, = 0, L„ = 2V,y, Ly = 0 

Using these expressions, it is not hard to see that 

n 

VF {v^,vy) = Y M^y"+^^yi' 

i=0 

□ 


C. Proof of lemma 

Proof Suppose 31U, Xi G 5^” such that the constraints are 
feasible for some p. Then, clearly, with the same choice, they 
are also feasible for any smaller p. Then, the sub-level sets 
of the problem are convex, and hence the problem is a quasi- 
convex optimization problem. Suppose XqLo 
is full-rank. When p = 0, for any fixed W, we can choose 
Xk = 0 and Ti,,Sk to satisfy the final equality constraints. 
Then, one needs to find W such that 


Sy W 


M,vf + N,vl 


\i=0 


F ml 


This can always be done if ^i'^f + Nivf is full rank, 
since one can simply choose 


lU = m > Mivf + N.yy 


Hence, this optimization problem is always feasible when 

z;r=o ^ 







